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Abstract. We present some numerical findings concerning the nature of the blowup 
£S) vs. global existence dichotomy for the focusing cubic nonlinear Klein- Gordon equa- 

tion in three dimensions for radial data. The context of this study is provided by the 
classical paper by Payne, Sattinger [32], as well as the recent work by K. Nakanishi, 
and the second author [29] . Specifically, we numerically investigate the boundary of 
£S) the forward scattering region. While the results of [29] guarantee that this boundary 

is smooth at energies which are near the ground state energy, it is currently unknown 
whether or not it continues to be a smooth manifold at higher energies. While we 
do not find convincing evidence of either smoothness or singularity formation, our 
numerical work does indicate that at larger energies the boundary becomes much 
43 more complicated than at energies near that of the ground state. 

B 

1. Introduction 

> 

This paper is concerned with the nonlinear Klein-Gordon equation 

^ nu + u = u tt - Au + u = u 3 , (t,x)ER 1+3 (1.1) 

with radial data. This equation exhibits several conserved quantities, of which the 
O energy, with u = (u, ii), 

> m= [ \^{\Vu\* + u* + i?)-\\ut]dx (1.2) 

is the most important one for our purposes. Note that we are only considering the 
focusing case here. Let us briefly review some basic well-known facts concerning the 



well-posedness of (1.1) in the energy space (this does not require any radial assump- 
tion) . 
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Theorem 1.1. The NLKG equation (1.1) is locally wellposed for data in if 1 (IR 3 ) x 
L 2 (R 3 ) in the usual (Duhamel formulation) sense and the energy is conserved. Small 
data lead to global existence and scattering to zero, whereas data of negative energy 
lead to finite time blowup in either time direction. If an energy solution exists on 
[0, T*) with T* < oo maximal, then necessarily ||w||i,3([o,r»);Lg) = °°/ if T* = oo and 
||ii||i3Qo,cxj);L|) < °°; then u scatters to zero as t — > oo. Finally, smooth data lead to 
smooth solutions. 

The proof is a straightforward application of energy and Strichartz estimates, see 
for example |39j. In the defocusing case one has global existence and scattering to 
zero for all data, see [21], (Zj, [E], pi], [IE], [33]. 

In what follows, "H := iJ 1 (IR 3 ) x L 2 (IR 3 ) is the phase (energy) space, and || • ||st stands 



for the L 3 L^-Strichartz norm associated with (1.1), respectively, on the time-interval 



[0, oo). Let us now define the forward scattering set. 
Definition 1 . The forward scattering set is defined as 

S+ := j(u ,ui) G U | u(t) := S(t)(u ,ux) 3 Vt > 0, ||«|| ST < ooj (1.3) 

where S(t) is the nonlinear evolution of the NLKG equation. 

The following properties of S + follow from the main wellposedness theorem above 
and some simple perturbative arguments. 

Lemma 1.2. The set S + enjoys the following properties: 

• S + D B$(0), a small ball in % 

• S+ is an open set 

• S + is path- connected. 

For the third property one can use that solutions of negative energy blowup in finite 
time, see [28], [32]. The following questions pose themselves quite naturally now: 

(1) Is S + bounded in T-Ll 

(2) What does dS + look like, is it a smooth manifold, or can it be very rough? 

(3) If dS + is a smooth manifold, does it separate a region of finite time blowup 
from one of global existence, at least locally? 

(4) What is the dynamical behavior of solutions starting on dS + l Are there any 
special solutions with data on dS + l 



Some partial answers to these questions were found in [29 J for ( 1.1 ) with radial data 
(the dim = 3 case with nonradial data was treated in [31], the cubic NLS equation 
in dim = 3 in [30], and the critical wave equation in dim = 3, 5 in [25], whereas the 
one-dimensional case is studied in [26J). To be more specific, it was found that the 
answer to (1) is "no". In fact, there exists a curve in % connecting to oo and such 
that a solution starting from an arbitrary point on that curve belongs to S + . The 
answers to the other questions are formulated in terms of the ground state Q. This 
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refers to the unique positive, radial, stationary energy solution of the NLKG equation. 
In other words, Q = Q(r) > 0, Q G H 1 (R d ), and 

-AQ + Q = Q P (1.4) 



Amongst all solutions ip of ( |1.4[ ), the ground state minimizes the energy 

1 -\ip\ p+1 )dx 



J{fp) :-- 



(\i\^\ 2 + V 2 ] 



P + 



In dim = 1, (1.4) has only two solutions (up to translation symmetry) which decay 
at ±oo, namely 



Q(x)=acosh p(f3x), a 



,v + l,^— n p — 1 



;i.5) 



In contrast, in dim = 3, (1.4) has infinitely many H 1 solutions, but up to translation 
symmetry only one positive solution, namely the ground state, see [38] . [I], [H]. In 
their classical paper [32], Payne and Sattinger showed the following, with 



K(<p):=d x J{<*<p) = [\Vv\ 2 + W\ 2 -\v\ P+l ](x)dx 

A=0 J R d 

being the scaling functional. 

Theorem 1.3. The two regions 

PS+ := {(uo,Wi) e U | S(uo,ui) < J(Q), K{u ) > 0} 
PS_ := {(u , Ul ) e H | E(u , Ul ) < J(Q), K(u ) < 0} 



+ ; 



are invariant under the flow of (1.1) in the following sense: if (m(0),m(0)) G PS 
then (u(t), u(t)) G PS + for as long as the solution exists, and the same holds for PS_. 
Solutions of (1.1) which lie in PS + exist for all times, whereas those in PS_ blow up 
in finite time (in both temporal directions). In particular, data of negative energy blow 
up in finite time, and S + y^T-L. 

The invariance is a fairly immediate consequence of the fact that the variational 
problem 

M{J(ip) \ peH 1 \{0},K&) = 0} 
has ±Q as the unique solutions (up to translation). Furthermore, in PS+ one has 

E{u) ~ \\u\\ H 



whence Theorem 1.1 implies global existence. The blowup in PS_ follows from a 
concavity argument. Very recently, Ibrahim, Masmoudi, Nakanishi [20J proved that 
solutions in PS + scatter to zero as t — )■ ±oo. Thus, PS + C S + . Returning to describing 
the partial answers to questions (2), (3), (4) above, introduce 

n £ := {(«o, ui)en\ E(uo, u x ) < J(Q) + e 2 } 
Then the following was shown in [29], [26J: 

Theorem 1.4. For some small e > one has 

• dS + fl "H e is a smooth, co-dimension one manifold 
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• it does separate a region of scattering to zero as t — » oo from one exhibiting 
finite time blowup 

• any data from dS + PI T-i £ lead to global solutions in forward time which scatter 
to ±Q, i.e., 

u = ±{Q,0)+v + o H {l) t-^oo (1.6) 

where v is an energy solution of the free Klein- Gordon equation. 

These results are for radial data, see [31] for a version with nonradial data. Further- 
more, the aforementioned references contain a description of the boundary dS + R "H e 
as the center-stable manifold associated with (Q,0) in the sense of Bates, Jones pQ. 
Moreover, the stable (and unstable) manifolds are one- dimensional and are described 
in the terms of the threshold solutions found by Duyckaerts, Merle in the energy 
critical case 



A well-known consequence of Theorem 1.3 is the instability of Q. More precisely 



in any neighborhood of (Q, 0) in H we can find data which lead to finite time blowup 
as well as global existence, respectively. To see this, we remark that the linearized 
operator L + = —A + 1 — pQ p ~ x has a single negative eigenvalue, the ground state. It 
is evident that L + has negative spectrum since 

(L + Q\Q) < 0. 

Let L + p = —k 2 p where k ^ and ||p|| 2 = 1. For further discussion of the spectral 
properties of L + and L_ = —A + 1 — Q p_1 , see [21] for the one-dimensional case, 
as well as [291 Lemma 2.3], [13], [12] for the three-dimensional case. Expanding the 
stationary energy J and the functional K around Q yields 

AQ + v) = J(Q) + \(L + v\v) + 0(\\vf Hl ) 
K(Q + v) = -(p-l)(Q p \v}+0(\\v\\ 2 Hl ) 



as — > 0. Therefore, by Theorem 1.3, data (Q + ep, 0) blow up for e > small 



and scatter to zero for e < small. The same holds for (1 + e)Q. 

Expanding further, we set u = Q + \p + w where w _L p. Then u = Xp + w and 

E(u) = J{Q) + l -{\ 2 - k 2 \ 2 ) + ±((L+w\w) + \\w\\l) + 0(|M|^) (1.8) 

where v = Xp + w. Since (L + w\w) > 0, we conclude that the only way to decrease 
the energy below J(Q) is through —k 2 X 2 . Moreover, in a ^-neighborhood of (0, 0) in 
the (A, A)-plane the set {E{u) < J(Q)} looks like 

{^ v )eR 2 \e-v 2 <o, \t\ + \v\<6} (1.9) 

at least up to cubic corrections. 

Needless to say that there are other important aspects of the problem which we do 
not touch upon in the present work. For instance, it is known [19] that the blowup 
rate is determined by the ODE u tt = u 3 , very similar to the analogous situation for 
the wave equation 
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2. The numerical methods 

At the core of all numerical work lie second order finite difference schemes (both 
explicit and implicit schemes were used) which solve the cubic NLKG equation with 
radial data. Such schemes are quite standard and are discussed in more detail below. 



In order to study the set S + for (1.1) numerically, the authors applied these solvers to 
data of the form (Q + Af, Bg) or (Af, Bg) where /, g are fixed radial functions, such as 
Gaussians, exponentials, with or without oscillatory factors (typical ones are sin(6r) 
or sin(6r 2 )). Here (A,B) runs over a fine grid of points from a suitable rectangle. If 
numerical blowup is found, which simply means that the H 1 x L 2 norm on some fixed 
ball becomes very large in finite time, something like 10 5 times larger than the norm 
of the initial data, the pair (A, B) is stored in a blowup file; if the H 1 x L 2 -norm 
becomes very small on some ball (say, falls below 1/10 of the size of the data), and 
remains small, then the pair (A, B) is regarded as "dispersive". If no decision of this 
type can be reached up to some prescribed maximum number of time steps, the data 
set is characterized as "indecisive" . Note that the latter case occurs if (Q + Af, Bg) 
falls exactly on the center-stable manifold (an extremely unlikely event; however, one 
may come very near to this manifold, which then requires a large number of time-steps 
in order to reach the blowup/dispersion decision). 

While this characterization is somewhat delicate (for example, it is not apriori clear 
why solutions could not blow up on a sphere of large radius, cf. [31]), the authors have 
found empirically that it leads to consistent and meaningful results (more on this can 
be found in Section |3| . 

The numerical calculations were performed with difference schemes applied to the 
equation 

v u ~ v rr + v = ^v 3 (2.1) 

for v(t,r) = ru(t,r). The choice of the difference scheme is a somewhat delicate 
matter. For the free wave equation there are several schemes, accurate to second or 
higher orders, which are known to converge |41j . For the nonlinear case, especially 



with the singularity on the right-hand side of (2.1 ), there does not appear to be much 
known in terms of the convergence of difference schemes for (2.1). 

The first, implicit, second order accurate difference scheme which we use is due to 
Strauss, Vazquez [ID] (which they introduced for the defocusing equation) and reads 
as follows: 

vj +1 - 2v n 3 + v]~ l v] +1 - 1v n } + 1 



Wf 1 + vj- 1 } 



(At) 2 (Ar) 2 2' 

1 [(^ n+i ) 3 + (^vr 1 + vf 1 ^- 1 ) 2 + (v^- i f] = o 



(2.2) 



4(j'Ar) 



G«+ 1 )-G«- 1 ) 



the nonlinear term being — n _\ — with G(v) = —\v . Here, we use the standard 

shorthand notation := v(nAt,jAr). The discretization steps At, Ar are chosen 
such that ^ = 0.9; in particular, the Courant-Friedichs-Lewy condition is satisfied. 
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Moreover, typical values of Ar are on the order of 10 3 . For the implementation, 



one applies a Newton scheme to solve for v ™ +1 using as a starting guess the explicit 
version of (2.2) in which all but the first occurrences of vj +1 are replaced by t>™, say. 
Note that due to the focusing sign the Newton scheme may become degenerate, which 
means that the derivative of the polynomial of vj +1 as defined by (2.2) can vanish 
or almost vanish. However, one checks that this can only happen if M(t) ■ (Ar) > 1 
where 

M(t) = max \u(t, r)\ = max \v(t, r)|/r; 



r>0 



r>0 



if the latter case occurs, then we terminate with a decision of blowup. To calculate 
the solution of (2.1) on a rectangle (t,r) G (0,T) x (0,R) we use a space-time grid 
on the larger rectangle (0,T) x (0, R + T), since the latter is the smallest rectangle 
containing the domain of dependence of (0, T) x (0, R). 
The scheme (12. 21) conserves the discrete energy 



En 



Ar 



+ 



E 



At 



E 



Ar 



Ar 



E 



,n+l\2 



+ 



E 



+ (V 



n\4_ 



(2.3) 



8(jAr) 



and is stable and thus convergent |5T] for the free equation. 

The second, explicit, second order accurate difference scheme is as follows: 



v n + t 



2v 7 + v 'j 



n-l 



v? +1 - 2v] + v r ;_ 



,n\3 







, (2.4) 

(At) 2 (Ar) 2 3 (jAr) 2 v 1 ' K ' 

It is easier to implement, and runs faster since no Newton iteration is required. It 
is still stable and convergent for the free equation. While it does not conserve the 
discrete energy (2.3) exactly, the authors have found that it does so approximately in 
those examples which they considered. 

A theoretical comparison of the difference schemes (2.2) and (2.4) appears to be 
quite challenging. However, after running numerous computations with these schemes 
the authors have found that they lead to the same conclusions as far as the disper- 
sion/blowup dichotomy is concerned (for example, in terms of the pictures appearing 
in the following section). However, as far as the actual blowup solutions are con- 
cerned, significant differences can be found near the blowup time. But we reiterate 
that in all cases checked by the authors the qualitative behavior comes out the same, 
which means that the appearance of the sets in Sectio n [3| comes out the same. 

Other than running their computations with both (2.2) and (2.4) (the latter typi- 
cally with finer resolution due to its faster execution), and then repeating the calcu- 
lations with Ar/2, the authors have used the bisection method as a means of checking 
the validity of figures such as those appearing in Section [3| This refers to a process in 
which two data pairs d% := (Ai,B\) and c?2 := (^2,-82), the first leading to blowup, 
and the second to dispersion, are chosen in close proximity to each other. Then one 
bisects the line segment joining them, tests the midpoint for blowup/dispersion, and 
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then continues with the line segment connecting the midpoint to either d\ or cf 2 > so 
that the endpoints of this new segment exhibit opposite behaviors. Iterating this 
procedure until machine precision is reached, the authors have found that one obtains 
data pairs which "almost" lie on the center-stable manifold. This means that the time 
evolution of these data exhibit "metastable" behavior, or "quasinormal" ringing. The 
latter term is used in the theoretical physics literature which concerns itself with the 
decay of dispersive tails, as they appear for example in the context of the so-called 
Price law. See also Bizoh et al. [3], [BJ. 

Finally, as an additional check, the authors have also run computations with a code 
that uses a compactified radial coordinate. This approach is quite common in the 
physics literature and it allows one (at least in principle) to follow the evolution up 
to spatial infinity which is important to rule out the possibility of late time blowup. 
Fortunately, such an unpleasant phenomenon does not seem to occur and the Klein- 
Gordon equation under consideration is well-behaved in this respect. 

As far as the coding is concerned, the authors used the "C"-language. The compu- 
tations leading to the results of the following section were carried out on two machines: 
(1) a MacPro capable of running 16 processes in parallel. In that case, the authors 
relied on the "forkO "-command to distribute the (A, 5)-loop over 50 (typically) par- 
allel processes at any given moment. (2) The Kraken supercomputer (a Cray XT5) 
at NICS, which is part of the TeraGrid. 



3. The numerical findings 

In order to explore questions (2)- (4) above numerically, the authors investigated 
solutions generated by data sets 

Vi(f,9) :={(Q + Af,Bg) | \A\ < a, \B\ < b} 
V 2 (f,g) :={(Af,Bg) | |A| < a, \B\ < b} 

where (/, g) are some simple fixed radial rapidly decaying functions, such as Gaussians, 
exponentials, possibly multiplied with oscillatory factors, and a, b > are chosen on a 
case to case basis. In the examples shown below, /, g are smooth, including near the 
origin (in fact, they are analytic). Once a choice of (f,g) is made, the numbers, a, b 
are then determined by trial and error (weighing numerical cost against the desire to 
obtain a large enough section of S + and dS + ). 



Discretizing (A,B) by a sufficiently fine rectangular grid, and solving (1.1) at each 
grid-point by the difference schemes of the previous section, one then colors each grid- 
point according to whether the numerical solution appears to exhibit blowup (it is 
kept blank) or scattering to zero (colored red) for large time^j 

In this fashion, one obtains two-dimensional sections of S + . In addition, the authors 
super-imposed the Payne, Sattinger regions PS± over the resulting images so as to 



The figures shown below represent only a small portion of all data which the authors have 
produced. The selection shown here is intended to capture the most important features which can 
be seen from those data. 
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Figure 1. Numerical results for V 1 (e~ r2 , e~ r2 ) and V 2 (Q, Q) 



obtain a meaningful comparison to Theorem 1.3 The image on the left of Figure [T] 
shows the outcome of such a computation for the data choice 

(it(0), «(0))(r) = (Q(r) + Ae- r \Be- r2 ) 

with the horizontal axis being A, and the vertical being B. The green and blue regions 
are PS + and PS_, respectively. Note that they meet at (0,0), which is precisely the 
soliton Q, in the conic fashion described by (1.9). The red region are data that lead 



to global existence (it includes green as a subset), whereas those that appear white 
(which include blue as a subset) correspond to finite time blowup. The image on 
the right-hand side of Figure [l] is produced by data (AQ(r), BQ(r)). Figure [I] agrees 
with the results of |32j and [29]. Indeed, PS+ does appear as a subset of the red 
region, PS_ as a subset of the white set, and the boundary of the red region near 
(Q, 0) looks like a smooth curve. It continues to look like a smooth curve even far 
away from (Q, 0), but this is currently not backed by theory. While the computation 
of the green and blue Payne-Sattinger regions does not pose a serious obstacle (the 
energy J(Q) can by obtained once Q is determined by a shooting-method, say, using 
Mathematica), determining the red regions is of course the main difficulty For that we 
use the difference schemes described in Section |2j Figure [2] and [3] show the outcome of 
computations where the data do not pass through either (Q, 0) or (— Q, 0). A natura 
section of this type is T> 2 (e~ r2 , e~ r2 ), and it appears on the left-hand side of Figure [2 
The right-hand image of that figure has very similar features. In particular, "spikes" 
and a smooth boundary (we discuss these spikes in more details below). Note that 
in these figures we only consider positive values of A, which suffices by the reflection 
symmetry about the origin. The blue PS_ region does not intersect the rectangle 
chosen in the right-hand image of Figure [2j and thus does not appear in there. 

A common feature of Figures [TJ [2] is that the boundary of the red region appears 
as a union of smooth curves. A new feature in the form of "bubbles" is shown in 
Figure [3j At least judging from the resolution used in these images, they do not seem 



NUMERICAL STUDY OF NLKG 9 




to involve any singularities. We use the notation (r) = a/1 + r 2 in the definition of 
the data (we prefer (r) to r as the former is smooth in IR 3 , whereas the latter is not). 




1 2 3 4 5 6 1 05 1 1JS 2 2J 3 33 4 4.5 



Figure 3. Numerical results for V 2 {e w ,sin(6r 2 )e (r2 1 ) 2 / 10 ) and 
V 2 {r 2 cos(10r 2 )e-< r2 \ sin(6r 2 )e- 4 ( r2 - 1 ) 2 ) 

In contrast to these examples, Figures [4] and [5] display a rather different phenome- 
non, namely more and more bands emanating from the red (scattering) region as one 
zooms into small and smaller scales. To be more specific, the image on the right-hand 
side of Figure [4] as well as those in Figure [5] show three zooms with successively higher 
resolutions near the round tip of the left-hand image of Figure |4} 

With some patience the reader will not only be able to find the filaments that are 
visible in one scale again at the next finer scale, but also between those and the solid 
red region a very thin new filament emerging. Note that the green color in the three 
zooms signifies the "indecisive" points where the algorithm was not able to decide 
between blowup and global existence. It is this indecisive region that gets resolved at 
the next finer scale, with a new filament emerging from it. 
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I 2 3 4 S 3.1 3.2 3.3 3.4 3.5 3.6 3.7 



Figure 4. Numerical results for T> 2 (e < r > , e r ) and one zoom near the boundary 

The zooms were computed on the Kraken supercomputer at the NICS, which is 
part of the TeraGrid. The coarse appearance of the third zoom, shown on the right- 
hand side of Figure [4] is due to the limitations on the computing time available to the 
authors. 

At this point, it is completely unclear whether or not infinitely many of these 
"bands" emerge out of the scattering region. Of course, this would imply that dS + is 
not a smooth manifold — perhaps a somewhat counter-intuitive assertion. 

A very visible common feature of all figures of this section (excluding of course the 
zooms) are the thin "spikes" that emanate from the bulk of the red region (the forward 
scattering set dS + ). Even though current theory cannot touch such phenomena as 
these "spikes", we now offer some — admittedly daring — explanations as to why 
they might appear. First, let us assume the scattering hypothesis which states that 
solutions starting from data that fall exactly onto the boundary exist globally in 
forward time and scatter to either Q or —Q. In addition to the proof of this hypothesis 
in [29] for energies which are at most slightly larger than that of the ground state, 
circumstantial evidence gathered by means of the bisection method described in the 
previous section lends some credence to this conjecture. However, if there are indeed 
infinitely many filaments in Figures [3j [4] then it is very difficult to guess what the 
dynamical implications might be; in particular, it is hard to conjecture (or simulate 
for that matter) what kind of asymptotic states might be associated with such a 
complicated boundary. 

Be it as it may, let us assume the "scattering hypothesis". Then the boundary 
splits into two open sets which scatter to either Q or — Q under the nonlinear flow 
(for convenience, let us call these positive and negative boundaries, respectively). The 
openness of these sets is another leap of faith, but a natural one: it simply expresses 
that the scattering properties should be stable relative to small perturbations inside 
the boundary. If this is so then these sets will not meet up smoothly but rather avoid 
each other. If a two-dimensional section does indeed display smooth boundary curves 
as in the numerical pictures, then the positive and negative curves cannot approach 
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each other within a compact set. Therefore, they have to move off to infinity. Note 
that this does not exclude that two boundary curves form a "tube" of some finite 
thickness, rather than a spike. Such a tube, however, is extremely unlikely since 
it would insinuate something like a uniform stability property at larger and larger 
energies. 




Figure 5. Two further zooms near the boundary of Fig [4] 

In this context, it would be very instructive to consider the Schrodinger equation 
since the soliton manifold in that case is connected (whereas for the NLKG equation 
it is discrete: {±Q}). If the previous explanation has any merit, then one might 
expect that spikes do not form in the Schrodinger case. Some preliminary numerical 
experiments by the authors appear to confirm that. More precisely, two-dimensional 
sections of the scattering region for the cubic radial focusing NLS equation in IR 3 
appear as round and smooth shapes — at least in the very limited number of examples 
the authors computed for NLS. 

We remark that the "spikes" should not be confused with the one-dimensional 
stable and unstable manifolds emanating from ±Q (or rather, neighborhoods thereof). 
Indeed, those manifolds would be only visible in very carefully chosen curved two- 
dimensional sections passing through ±Q, and finding such a section would be a zero 
probability event in any reasonable sense. But this contradicts the fact that the spikes 
appear in each of the numerically computed images (but not in the zooms, of course) 
that the authors produced. 

Needless to say, it is essentially impossible to hit the boundary exactly just by 
sampling a grid of data as we do to produce the figures of this section. Nevertheless, 
one can test the aforementioned scattering hypothesis by means of the "bisection 
method" described in Section |2j see also Bizoh et al. [5] , [6] . This refers to taking 
two data pairs, one from the red region, and another from the white region which 
are in close proximity of each other. Then one numerically computes the solution 
starting from the midpoint between these two points, tests for blowup/dispersion, 
and continues bisecting the interval whose endpoints exhibit opposite behaviors. In 
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this way one obtains data that have a clear tendency of approaching either Q or 
— Q for large times, albeit with possible oscillations about these two states, followed 
by eventual blowup or dispersion. The aforementioned oscillations are referred to as 
"quasinormal ringing" in the physics literature, see for example Bizoh et al. 

To the authors, this bisection process also provides a strong check of the accuracy 
and validity of the numerically computed region S + . Indeed, if the initial data pair 
for the bisection process does not lie close to, and on opposite sides of, the actual dS + 
- which is what the figures claim — then the bisection process would not converge 
in the described fashion. 




Figure 6. P2(^sin(6r)e _ 2\ r / ) cos(6r)e~ (r ) and the curved section from (3.1) 



Figure [6] shows two more sections with analytic, decaying data. The left-hand image 
is a planar section as defined above, and it provides (after the reflection about the 
origin) an example with six spikes rather than two as in the previous examples. A 
priori, there seems to be no reason to assume that the number of spikes could not be 
arbitrarily large. 

The image on the right-hand side of Figure [5] is different from all those previously 
considered in so far as it is not planar. More specifically, it shows a section produced 
by 

(/, g) = {{A + B)B cos(2(A - B)r)exp(--^), (A - E)exp(-r 2 )) , (3.1) 

with \A\ < 6, \B\ < 4. 

The final series of numerical experiments which we present here concerns three- 
parameter families. This refers to data of the form (Af(C,r),Bg(C,r)), say, where 
A,B,C are parameters which are chosen in some finite range. For fixed C these 
data generate a two-dimensional section as shown in the previous images, but as C 
varies one obtains a whole family of sections which move "continuously" (at least, if 
C changes only by very little). The specific data leading to Figures [7j [8] are 

P 2 (r 2 cos(L7r)e- r2 ,r 3 sin(6r)exp(-2(r 2 - l) 2 ) (3.2) 
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Figure 7. A series of sections as in (3.2) with increasing C 



where C is in the vicinity of 9.72, 9.73 and up. Figure [7] shows a series of such two- 
dimensional slices obtained by letting C increase through six values (the white vertical 
lines visible in the red region are an artifact caused by a slightly coarser resolution 
in the horizontal direction). What can be seen in these six images, starting with the 
upper left, is first the merging of two disjoint sets, followed by the pulling apart of 
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two sets which then become disjoint. Moreover, the direction in which the sets pull 
apart is roughly orthogonal to that in which the merging takes place. 








Figure 8. A series of higher-resolution images for (3.2) with increasing C 



The images in Figure [7] were computed on the MacPro and use a fairly low resolution 
in order to keep the computation times reasonable. A much finer resolution was 
obtained on the supercomputer Kraken at the NICS, leading to the sequence of images 
shown in Figure [8j They are zoom-ins into the area where the pinching takes place as 
shown in Figure [7j and depict the process of merging and pulling apart more clearly. 
The uneven appearance of some of the upper boundaries in the images of Figures [8] is 
due to the fact that in most computations the so-called wall-time was reached which 
then leads to termination. 

Interpreting these images is a delicate issue, as they are two-dimensional sections 
of an infinite-dimensional object, namely S + . Note carefully that they certainly do 
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not imply singularity formation of the actual boundary dS + as can be seen from the 
example of the one-sheeted hyperboloid in K 3 , say with rotation axis given by the 
a^-coordinate direction. Then as the planes x\ = const slice through this hyperboloid 
they will produce an image given by x\ — x\ < at the moment of tangency to 
the waist line. This is what seems to happen between the fourth and fifth image of 
Figure |8j 

4. Concluding remarks 

The authors hope that this note will encourage further numerical, or analytical, 
investigations of the set S+, both for the NLKG equation in various dimensions and 
with various powers in the nonlinearity, as well as other nonlinear dispersive equations 
such as the nonlinear Schrodinger equation. It would be highly desirable to employ 
different, and possibly more sophisticated, numerical algorithms. For example, one 
could try the symplectic difference schemes developed by Duncan [H] and Reich [35] . 
or a different approach altogether. 

The most important problem arising in this context is to decide whether or not 
the boundary of the forward scattering region can indeed become singular at larger 
energies. And if so, in what ways the boundary may degenerate. As discussed in the 
previous section, the only evidence we have found that points towards the possibility 
of some sort of singular behavior is the emergence of more and more filaments from 
the scattering region. The next step would then be to see if any possible degeneration 
of the boundary is reflected by anomalies in the actual wave dynamics. 

Upon request the authors will email all numerical data which they have produced 
in the course of this work (which is considerably larger than that presented here). 
However, hopefully other researchers will carry out their experiments with data and 
numerical solvers of their own choosing. It seems most important to conduct a sys- 
tematic and exhaustive study of the boundary of <S+ by means of the bisection method 
which we discussed above in order to test the "scattering hypothesis" . In the process 
one can also test the accuracy of the numerically computed two-dimensional sections, 
particularly if a different solver is employed for the bisection than the one used for 
the computation of the two-dimensional slices of S + . 

More work needs to be done on three-parameter families. In particular, it would 
be of interest to study how the "bubbles" in Figure [3] evolve as a third parameter is 
being changed. Furthermore, one could use such a three parameter family to obtain a 
two-dimensional section constant energy. And thirdly, it would be very appealing 
to produce three-dimensional images of the intersections of <S+ with a cube in R 3 (this 
requires substantial computational resources). 

The issue here is not so much to produce more pictures, but to pinpoint various 
phenomena relating to the topology of dS + and then to connect them to the actual 
wave dynamics. An example of this is the aforementioned scattering hypothesis but 
one can also ask about the dynamical implications (if any) of the filaments as they 
appear in Figures [4] and |5j 



'We thank Kenji Nakanishi for this suggestion. 
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Finally, it would be interesting — as well as challenging — to consider the same 
questions that we pursue here for other equations, such as the one-dimensional NLKG, 
the Schrodinger, as well as the energy critical wave equations (i.e., with u 5 nonlin- 
earity) in IR 3 . The one-dimensional equation has considerably less dispersion than 
the three-dimensional one, so it is perhaps somewhat more delicate to rely on the 
numerics in that case. In this context, one can also investigate the dependence on the 
power of the nonlinearity, see Bizon [6] for such considerations in connection with the 
aforementioned "quasi-periodic ringing" and the emergence of a threshold resonance 
at power p = 3. It is natural to consider even data as in but general data are 
also important, of course, and considerably more challenging numerically. As for the 
Schrodinger and critical wave equations, these are numerically (as well as analytically) 
harder and the authors have done only preliminary exploratory work on them as far 
as the numerics is concerned. 
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